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Abstract 

A numerical simulation of the incompressible viscous flow through a prosthetic 
tilting disk heart valve is presented in order to demonstrate the current capability 
to model unsteady flows with moving boundaries. Both steady-state and unsteady 
flow calculations are done by solving the incompressible Navier-Stokes equations in 
three-dimensional generalized curvilinear coordinates. In order to handle the moving 
boundary problems, the chimera grid embedding scheme which decomposes a complex 
computational domain into several simple subdo ma i ns is used. An algebraic turbulence 
model for internal flows is incorporated to reach the physiological values of Reynolds 
number. Good agreement is obtained between the numerical results and experimental 
measurements. It is found that the tilting disk valve causes large regions of separated 
flow, and regions of high shear. 

Introduction 

Various types of prosthetic heart valves have been used widely as the replace- 
ments of natural valves since the first successful valve replacement performed in 1960. 
However, each of the valve design has some difficulties, which cause the artificial heart 
valve to be less efficient than the natural one. The difficulties which are related to 
the nonphysiological flow characteristics of the currently used prosthetic heart valves 
are: 1) Large pressure losses across the valves prevent the heart working efficiently; 2) 
Separated and secondary flow regions cause clotting; 3) High turbulent shear stress can 
damage the red blood cells. Having detailed knowledge of the flow quantities can help 
a design engineer improve the valve geometry, where a smooth flow is desired. Certain 
experimental studies 1 _s have pointed out the effects of the stagnation and recircula- 
tion regions and compared commonly used valve geometries. Since the experimental 
measurements provide flow characteristics for only certain regions of the flowfield, the 
numerical simulation of the flow through the artificial heart valve will be extremely 
helpful in the design and development stage of the prostheses. 

Most of the numerical studies modeled the flow through the heart valve devices 
by excluding the moving boundary problems. Underwood and Mueller 4 obtained the 
flow characteristics for Kay-Shiley disk type valve using the stream function-vorticity 
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formulation. Their results showed agreement with experimental data up to Reynolds 
number of 600. Idelsohn, Costa, and Ponso 5 modeled the flow through Kay-Shiley 
caged disk, Starr-Edwars caged ball, and Bjork-Shiley tilting disk valves and com- 
pared their performance. A maximum Reynolds number of 1500 was reached in their 
numerical study. In the above studies, the caged disk and caged ball geometries are 
axisymetric, and the tilting disk geometry is two-dimensional. In actual case, the tilt- 
ing disk geometry is three-dimensional, the flow through heart valves is unsteady, and 
Reynolds numbers are as high as 6500. Peskin and McQueen 8 modeled the prosthetic 
heart valves in the numerical simulation of the flow in the heart. They used boundary 
forces derived from the energy function in order to model valve opening and closing, 
and they also modeled the elastic behavior of the walls. Their solution is obtained for 
the low Reynolds numbers in two dimensions using square cartesian mesh. McCracken 
and Peskin 7 applied a combination vortex-grid method for the blood flow through the 
mitral valve in two dimensions. This method is applied to the problems in which 
the solution does not have strong dependence on the Reynolds number. Peskin and 
McQueen 8 demonstrates the capability of modeling elastic behavior of the heart muscle 
by applying their extended three dimensional solution procedure to a toroidal tube. In 
order to obtain a solution procedure turned at design improvements in prosthetic heart 
valves, the computation of steady-state and unsteady flow through the Bjork-Shiley 
tilting disk valve in three-dimensional configuration with the use of a grid embedding 
scheme is proposed in the current work. The equations are solved a in curvilinear 
generalized coordinate system, and the valve opening and closing are simulated by 
calculating the forces acting on the valve. 

One of the biggest difficulties in the simulation of the flows in complicated 
three-dimensional configurations is the discretization of the physical domain with a 
single grid. The problem becomes more severe if one body in the domain of interest 
moves relative to another one as is seen in the tilting disk configuration. The use of 
a zonal approach 9 would be a practical solution of the moving boundary problem if 
the grids could be constrained to common boundaries. The chimera grid embedding 
technique 10 ’ 11 provides a greater flexibility for the grid motion. Instead of using com- 
mon boundaries between grids, common regions are used. In the present work, the 
chimera approach is used to discretize the geometry of the disk valve. In addition, the 
procedure obtained for the heart valve configuration can be easily utilized for other un- 
steady incompressible viscous flows with moving boundaries, e.g., flow through Space 
Shuttle ExternaFIank/Orbiter propellant feed line disconnect flapper valves. 

In the first section, the method of solving the Incompressible Navier- Stokes equa- 
tions is described, and the algebraic turbulence model is summarized. Next, the geom- 
etry ancTthe use of the chimera scheme are discussed. Following that is a presentation 
of the computed results obtained from the current approach*— ~ 

Governing Equations and Method of Solution 

The algorithm used in both steady-state and unsteady flow calculations is based 
on the method of artificial compressibility, which produces a hyperbolic system of 
equations by introducing a time derivative pressure term into the continuity equation. 
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The resulting incompressible Navi er- Stokes equations can be written in a generalized 
curvilinear coordinate system (£,*7?0 follows 

* +T( {E - e ’ )+ & f - f ’ )+ tS g - g ’ )=0 (1) 

where Q, and convective flux vectors E, F , G are 
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Here J, P, p, u, r, and w denote the Jacobian of transformation, the pseudocompress- 
ibility coefficient, pressure, and velocity components, respectively. The contravariant 
velocity components U , V , and W are defined as 
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For an orthogonal grid assumption, the viscous flux vectors E v , F v , and G v are given 
by 
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where Rt is the Reynolds number. 

In the steady-state formulation, the time derivatives are differenced using the Eu- 
ler backward formula. The equations are solved iteratively in pseudo-time until the the 
solution converges to a steady state. In the time-accurate formulation, the time deriva- 
tives in the momentum equations are differenced using a second-order, three-point, 
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backward-difference formula. The equations are iterated to convergence in pseudo- 
time for each physical time step until a divergence-free velocity field is obtained. The 
numerical method uses a second-order central difference for viscous terms and a higher 
order flux-difference splitting for the convective terms. The £ derivative of the convec- 
tive flux E can be written as 


dE _ [E{+ 1/2 - #1-1/2] 
dt ~ A* 

The numerical flux Ei+i/ 2 defined as follow 

-E*+ 1/2 = - [-®W»+i ) + E(Qi ) — &+1/2] ( 2 ) 

where the ^t+1/2 is a dissipation term. The order of the scheme is determined by the 
definition of the dissipation term &+1/2* For &+1/2 = the differencing is reduced to 
a second-order central difference scheme. A first-order upwind flux is defined by 



and a third-order upwind flux is given by 

&+1/2 = ~ 1/2 — ^^H-i/2 + ^^t+1/2 “ ^^i+ 3 / 2 ) W 

where AE* is the flux difference across positive or negative traveling waves, and is 
computed as 

here A* is the plus (minus) Jacobian matrix. The A operator, and Q are given by 

AQi+1/2 = Qi+i — Qi 

Q=\{Qnrl+Qi) 

An implicit delta law form approximation to Eq.(l) after linearization in time 
and the use of approximate Jacobians of the flux differences results in a 4 x 4 block 
heptadiagonal matrix equation. The matrix equation is solved iteratively by using a 
nonfactored line relaxation scheme, which maintains stability and allows a large pseudo- 
time step to be taken. At each sweep direction, a tridiagonal matrix is formed and off 
line terms of the matrix equation are moved on the right-hand side of the equation. 
Details of the numerical method are given in Refs. 12 - 14 . 

An algebraic mixing-length turbulence model, which is presented in Ref. 15 , is 
utilized in the present computations. The turbulent eddy viscosity is taken as 

v t = / J H 
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where l is the mixing-length, and |u>| is the absolute value of vorticity. In order to 
account the effect of more than one wall in the region close to the tilting disk, the 
mixing-length is given as 

H=1 

where n is the number of walls, c = 0.4 is the Von Karman’s constant, yu is the distance 
from the k th surface, and Dk is the Van Driest damping factor for the k th surface. 

Geometry and Grid System 

In the Bjork-Shiley tilting disk heart valve, the tilting disk is placed in front 
of the sinus region of the human aorta. The aortic root has three sinuses about 120 
degree apart from one another. The tilting disk valve model used in this computation 
is simplified by assuming that the sinus region of aorta has a circular cross-section. The 
cage and struts which hold the free-floating disk inside the sewing ring are not included 
in the geometry. It is also assumed that the walls do not have an elastic deformation. 
The computational geometry used in unsteady flow computations is given in Figure 
1. The channel length is taken to be five aorta diameters long. The disk motion is 
illustrated by showing three different positions of the disk. The disk angles shown are 
75, 50, and 30 degrees as measured from the centerline of the aorta. The tilting disk 
is allowed to rotate about the horizontal axis that is 1/6 of a disk diameter below 
the center of the disk. Because of this asymmetric disk orientation, the flow is three 
dimensional. 

The chimera grid embedding te chni que, which has been successfully used for 
external flow problems, has been employed by using two overlapped grids as shown in 
Figure 2-a. Grid 1 contains 17,199 points which are distributed 63 x 21 x 13 in the (, q, 
and (, directions. It occupies the whole region in the aorta from entrance to exit, and 
remains stationary. Grid 2 has 4,725 points as distributed 25 x 21 x 9. It wraps around 
the tilting disk, and moves with the disk. The lateral symmetry planes of the two grids 
are shown in Figure 2-b in order to demonstrate how the grid embedding scheme is 
applied to the present problem. Grid points which lie within the disk geometry and 
outside the aorta grid are excluded from the solution process. These, excluded points 
are called hole points, and the immediate neighbors of the hole points are called fringe 
points. The information is passed from one grid to another one via fringe and grid 
boundary points by interpolating the dependent variables. Tri-linear interpolation is 
used in the present computations. In order fo distinguish the hole and fringe points- 
from regular computational points, an IBLANK array is used in the flow solver. For 
hole, grid boundary, and fringe points IBLANK is set to zero, otherwise it is set to one. 
In order to exclude the hole and grid boundary points from the solution procedure at 
each time step, the A Q solution is processed as follows 

Qn+i _ q» + * IBLANK 

Since 3”* order flux-differencing is used for convective terms, the order of differ- 
encing needs to be reduced to the second order differencing near the fringe points. 
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Computed Results 


Presented here are the results of steady-state flow and unsteady flow with 
the disk motion in the configuration mentioned above. The problems are non- 
dimensionalized by using the entrance diameter as a unit length, and the average inflow 
velocity as the unit velocity. The geometries used in the steady-state and unsteady cal- 
culations are similar to the geometries used in the experimental studies in Ref. 1,2 
and Ref. 3 respectively. In order to reduce the computational effort and memory size, 
the inflow and outflow boundaries are placed a short distance from the region of in- 
terest in comparison with the boundaries in the experimental studies. In addition, the 
exact shape of the sinus region of aorta used in the experiments is not available at 
present. Because of these differences between experimental and computational config- 
urations, there are small differences between experimental measurements and present 
computations. 

At a viscous no-slip surface, both the velocity and the pressure gradient normal 
to the wall are specified to be zero. At the inflow boundary, the velocities are specified, 
and the pressure is determined from the characteristic boundary condition. At the 
outflow boundary, static pressure is specified, and velocities are calculated from the 
three characteristic waves traveling out of the computational domain. 

Steady-state calculations for the 30 degree disk orientation have been carried 
out for Reynolds numbers in the range of 2000 to 6000, in which experimental data is 
available. The Reynolds number is based on the diameter and the mean velocity at 
the entrance of the channel. Figure 3 shows the convergence history for a Reynolds 
number of 5972. Both averaged residual and maximum divergence of velocity have 
dropped ten order of magnitude in 600 iterations. For an overlapped grid application, 
the convergence is shown to be very fast. The values for grid 1 are drawn with solid 
lines, and the values for grid 2 are drawn with dashed lines. The computing time per 
grid point per iteration is about 2xl0~ 4 sec. Figure 4 shows the pressure drop across 
the Bjork-Shiley tilting disk valve at different flow rates of physiological interest. PI 
and P2 are the pressures at the points located 150 mm, and 20 mm upstream from 
the disk at the centerline of the channel respectively. In order to compare numerical 
results with the experimental measurements given in Ref. 2, the numerical results are 
redimensionalized. The computed and measured axial velocity profiles at 42 mm down- 
stream from the disk are shown in Figure 5. Figure 5-a demonstrates the horizontal 
plane, where axial velocity profiles are given, through the center of the channel. Axial 
velocity profiles given in Figures 5-b through 5-d illustrate how the stagnation region, 
which is created by the tilting disk valve, is dominant in the region of 1.5 disk diam- 
eter downstream. The numerical results are shown with dots and the experimental 
results are shown with triangles. The numerical results compare favorably with the 
experimental measurements given in Refs. 1-2. 

Velocity vectors on the lateral symmetry plane are given in Figure 6 for a 
Reynolds number of 5972. The flow, which is directed to the upper part of the aorta, 
generates vortices in the sinus region of the aorta and a large separated region along 
the lower wall of the aorta. Since separated and low flow regions have potential for 
thrombus formation, clotting may occur on the upper sinus region and the lower wall of 
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the aorta. The flow is highly accelerated new the tilting disk and the upper wall. Fig- 
ure 7 shows vortiaty magnitude contours on the surface of the channel, inflow surface, 
and outflow surface of the disk, respectively. It is assumed that maximum vorticity 
magnitudes indicate the regions of high shear. The sewing ring surface and the edges 
of the disk are the regions having maximum vorticity magnitude. The upper wall of 
the channel also has considerably high vorticity magnitudes. 

Unsteady flow calculations have been carried out in order to demonstrate and 
analyze the flow during disk opening and closing. For the present computation, one 
cycle of valve opening and closing requires 70 physical time steps. During each time 
step, subiterations are carried out until maximum divergence of velocity and maximum 
residual drop below 10 -s . The computing time required for one cycle of the valve 
opening and closing is aproximately 5 hr. During the valve opening, inflow velocity 
is imposed at the entrance of the channel. The inflow velocity is chosen as a sinuous 
function in time. The forces acting on the disk are calculated, and the disk rotation 
angle is determined. For large disk rotation angle, some information may be lost 
between the grids when the grid embedding technique is used. In order to prevent the 
information loss, the maximum allowed disk rotation angle at each physical time step 
is taken to be less than three degrees. As soon as the disk reaches its fully opened 
position, which is 30 degrees measured from the horizontal plane, the flow direction is 
reversed by imposing the inflow velocity at the exit of the channel. 

Figures 9-a through 9-f illustrate the velocity vectors on the lateral symmetry 
plane at t/T = 0.13, 0.285, 0.385, 0.53, 0.685, and 0.8 respectively. T is a period of 
one cycle during the valve opening and closing. The velocities are very high in the 
region between the disk and the channel wall as shown in Figure 9-a. During the disk 
opening, two vortices axe formed at the upper and lower edges of the disk. The flow 
starts to separate behind the disk and reattaches to the wall as shown in Figure 9-b. 
The stagnation region behind the disk moves downstream as the disk rotates. Highly 
skewed velocity profiles are seen downstream from the disk as illustrated in Figure 9-c. 
The growth of the vortices has also been observed m the sinus region of the aorta while 
the flow opens the valve. Along the lower wall a separation region is formed. Figure 
9-d shows the beginning of the valve closing. At this moment, the location of imposing 
the inflow velocity is moved from entrance to exit. Major flow near the upper wall of 
the channel forms a recirculation region downstream from the disk. With the help of 
this recirculation, the lower wall of the channel becomes the major flow region during 
the valve dosing,, and upper wall region becomes the low flow region. 

The vortiaty magnitude contours on the surface of the channd at t/T = 0.13, 
0.285, and 0.385 are shown in Figures 8-a through 8-c in order to indicate' the regions 
of high shear. At the beginning of valve opening, jet-like flow between the sewing ring 
and the disk cause very high vortiaty magnitudes as shown in Figure 8-a. During the 
disk rotation, the high vortiaty region on the upper wall of the channel moves from 
the sinus region to downstream as seen in Figures 8-a through 8-c. The results of more 
realistic flow calculations, such as the flow through the Pennsylvania State artifidal 
heart including tilting disk valve opening and dosing, will be reported in the future. 
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Concluding Remarks 


The solution procedure for unsteady incompressible viscous flow computations 
has been extended with the incorporation of the grid embedding approach. This has 
been used to simulate both steady-state and unsteady flow through a tilting disk valve. 
The physiological values of the Reynolds number have been achieved with the use of 
a simple mixing-length turbulence model. The numerical results for 30 degree dis k 
orientation were compared against the experimental data, and good agreement was 
obtained. The flow during the disk opening and closing were simulated within a rea- 
sonable computing time. The present capability of simulating complicated internal 
flow problems with moving boundaries is demonstrated. The procedure obtained here 
is quite general and applicable for various types of complicated geometries. 
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Figure 1. Tilting disk geometry showing valve opening. 




Figure 2. b) Side views of two overlapped grids before and after the hole points are excluded. 
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Figure 3. Convergence history for Re — 5972. 
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b) Re = 2390, Steady flow rate = 167 cm 3 /sec. 





ioa ia* 2QO 23 


DISTANCE FROM INSIDE WALL (mm) 

c) Re * 3380, Steady flow rate = 250 cm 2 /sec 


StMdy Flow Rate. Q (emOWe) 


Figure 4. Pressure drop across the tilting disk 
valve versus steady*state flow rate. 
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d) Re = 5972 , Steady flow rate = 417 cm 3 /see. 

Figure 5. Top view of horizontal plane and axial velocity profiles 
on the horizontal plane for different flow rates. 
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Figure 6. Side view of velocity vectors on the vertical plane. 
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Figure 7. Voriicity magnitude contours on the a) channel surface, b) inflow surface of the disk, 
c) outflow surface of the disk. 
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c) t/T = 0.385 

Figure 8. Vorticity magnitude contours on the channel 
surface during the valve openning. 


b) t/T = 0.257 
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f) t/T = 0.857 


plane showing the valve opening and closing 


